set more off
set maxiter 100
capture log close
clear
cd "/Users/danielhill/Documents/ZiOP_jcr"
log using "opsims4Jan13.log", replace

/* program that generates data w/ dv following ziopr dgp (called in loop) */

capture program drop syndata
program syndata, rclass
	args obs b0 b1 b2 g0 g1 tau_1 tau_2
	clear
	set obs `obs'
	gen x1 = -5+(10)*runiform()
	gen x2 = rnormal(0,2)
	gen z1 = rbinomial(5,.5)
    gen e1 = rnormal() 
    gen e2 = rnormal() 
	gen infl_lat = scalar(g0)+scalar(g1)*z1+e2
	gen infl = infl_lat > 0 
	sum infl, meanonly
	scalar _infl = r(mean)
	sum x1, meanonly
	scalar _meanx1 = r(mean)
	sum x2, meanonly
	scalar _meanx2 = r(mean)
	sum z1, meanonly
	scalar _meanz1 = r(mean)

	gen out_lat = scalar(b0)+scalar(b1)*x1+scalar(b2)*x2+e1
	gen outcome = out_lat > scalar(tau_1)
	replace outcome = 2 if out_lat > scalar(tau_2)

	gen y = infl*outcome
end

/* Create file and variables to store results from simulations */

postutil clear
postfile opsim obs sims b0 b1 b2 g0 g1 t2 /// 
		     ob0 ob1 ob2 ot2 ///
		     ovcv11 ovcv21 ovcv31 ovcv41 ///
		     ovcv22 ovcv32 ovcv42 ///
		     ovcv33 ovcv43 ///
		     ovcv44 ///
             meaninfl meanx1 meanx2 meanz1 rco ///
		     using "opsims4Jan13.dta", replace

/* Generate macros for number of observations in data, number of simulations, seed value*/

timer on 1
local obslist "4000"
local n = 5000
local sd = 666
local constant -6 -4.4 -3.3 -2.5 -1.7 -0.5 3 

/* Simulation loop */
	foreach obs of local obslist {
      local sd = `sd'
      qui set seed `sd'

      dis _n(2) as result "Number of Observations: `obs'"
      
	    foreach j of local constant {
      
			local i = 1

      		while `i' <= `n' {
	
			scalar b0 = 1
			scalar b1 = .5
			scalar b2 = 1.5
			scalar g1 = 1
			scalar g0 = `j'
			scalar tau_1 = 0
			scalar tau_2 = 2

		   	qui syndata `obs' scalar(b0) scalar(b1) scalar(b2) ///
        	scalar(g0) scalar(g1) /// 
			scalar(tau_1) scalar(tau_2)
	
			capture myoprob y x1 x2
			matrix eb = e(b)
			matrix ev = e(V)
			local rco = _rc
			if `rco'==0 {
			scalar ob0 = eb[1,3] 
			scalar ob1 = eb[1,1]
			scalar ob2 = eb[1,2]
			scalar ot2 = eb[1,4]
			scalar ovcv11 = ev[1,1]
			scalar ovcv21 = ev[2,1]
			scalar ovcv31 = ev[3,1]
			scalar ovcv41 = ev[4,1]
			scalar ovcv22 = ev[2,2]
			scalar ovcv32 = ev[3,2]
			scalar ovcv42 = ev[4,2]
			scalar ovcv33 = ev[3,3]
			scalar ovcv43 = ev[4,3]
			scalar ovcv44 = ev[4,4]
      		}
			else {
			scalar ob0 = . 
			scalar ob1 = .
			scalar ob2 = .
			scalar ot2 = .
			scalar ovcv11 = .
			scalar ovcv21 = .
			scalar ovcv31 = .
			scalar ovcv41 = .
			scalar ovcv22 = .
			scalar ovcv32 = .
			scalar ovcv42 = .
			scalar ovcv33 = .
			scalar ovcv43 = .
			scalar ovcv44 = .
      		}
	
/* Post results */

			post opsim (`obs') (`n') (scalar(b0)) (scalar(b1)) (scalar(b2)) ///
			(scalar(g0)) (scalar(g1)) (scalar(tau_2)) ///
            (scalar(ob0)) (scalar(ob1)) (scalar(ob2)) (scalar(ot2)) ///
			(scalar(ovcv11)) (scalar(ovcv21)) (scalar(ovcv31)) (scalar(ovcv41)) ///
			(scalar(ovcv22)) (scalar(ovcv32)) (scalar(ovcv42)) ///
			(scalar(ovcv33)) (scalar(ovcv43)) ///
			(scalar(ovcv44)) ///
            (scalar(_infl)) (scalar(_meanx1)) (scalar(_meanx2)) (scalar(_meanz1)) ///
            (`rco') 
            	
			drop x1 x2 z1 e1 e2 infl_lat infl out_lat outcome y
			if mod(`i',1000)==0 display as result "`i' " _c
 
			local i = `i' + 1
  			}
		}
	}
postclose opsim	
set more on
log close

timer off 1
timer list 1

